home *** CD-ROM | disk | FTP | other *** search
/ Gold Medal Software 3 / Gold Medal Software - Volume 3 (Gold Medal) (1994).iso / stats / lsqrft15.arj / FITCMDS2.C < prev    next >
Text File  |  1994-02-11  |  9KB  |  340 lines
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <string.h>
  4. #include <math.h>
  5. #include "fit.h"
  6.  
  7. extern int debug;
  8.  
  9. /* function to get data and return number of data points */
  10. int get_data(data, command, num_indep, inbuf, order,filename,maxrows, maxcols)
  11. double **data;
  12. char *command;     /* command which was typed at the fit> prompt */
  13. int num_indep;
  14. char *inbuf;           /* inbuf is really output buffer */
  15. struct data_order *order;
  16. char *filename;
  17. int maxrows;
  18. int maxcols;
  19. {
  20. int i, j, jmax, ndata;
  21. char str[20][30];
  22. FILE *instream;
  23. int iopt1;
  24. char  inbuf1[512];
  25. iopt1 = 20;
  26.  
  27. sscanf(command,"%s %s %d", inbuf, filename, &iopt1);
  28.  
  29. if(( instream=fopen(filename,"r") ) == NULL) return 0;
  30. i = 0;
  31. while( fgets(inbuf1,512,instream) != NULL && i < maxrows){
  32.     if(inbuf1[0] == '#') continue;
  33.     jmax=sscanf(inbuf1,"%s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s",
  34.             str[0],str[1],str[2],str[3],str[4],str[5],str[6],str[7],str[8],str[9],
  35.             str[10],str[11],str[12],str[13],str[14],str[15],str[16],str[17],str[18],str[19]);
  36.     for(j=0; j < jmax && j <= iopt1 && j < maxcols-3; j++){
  37.         data[j][i] = atof(str[j]);
  38.     }
  39.     i++;
  40. }
  41. fclose(instream);
  42. ndata = i;
  43.  
  44. /* assign first columns to be independent variables */
  45. for(i=0; i < num_indep; i++)order->x[i] = i;
  46.  
  47. /* assign next column to be dependent variables */
  48. order->y = num_indep;
  49.  
  50. /* if there is a column num_indep+1, assign it to be error in y */
  51. if(jmax > num_indep) order->isig = num_indep + 1;
  52.  
  53. /* assign a data column to be all ones: no weighting */
  54. for(i = 0; i < ndata; i++)    data[jmax][i] = 1;
  55. order->nsig = jmax+0;
  56.  
  57. /* assign a data column to be used for statistical weighting */
  58. order->ssig= jmax+2;
  59.  
  60. /* assign a data column to be used for yfit */
  61. order->yfit= jmax+1;
  62.  
  63. if(jmax == 2 && num_indep == 1) order->sig = 2;
  64. if(jmax == 3 && num_indep == 2) order->sig = 3;
  65.  
  66. /* Assigning data columns for different weightings */
  67. /* takes up extra memory, but is faster and simpler */
  68. /* especially in the case of statistical weighting: */
  69. /* we don't have to take the sqrt of y at every point */
  70. /* for every iteration.  It is done once when statistical */
  71. /* weighting is selected */
  72.  
  73. return ndata;
  74. }
  75.  
  76. /* This function initializes the parameters */
  77. void ip(command, inbuf, func, a, ma)
  78. char *command, *inbuf;
  79. int (*func)();
  80. double *a;
  81. int ma;
  82. {
  83. int i;
  84. char cmdargs[NUM_ARGS][30];
  85.  
  86. if(func == NULL)
  87.     printf("ip failed, you must select a function first\n");
  88. else if(parse(command,cmdargs) == ma){
  89.     for( i = 0; i < ma; i++){
  90.         /* a * in place of a number leaves the parameter unchanged */
  91.         if(strcmp(cmdargs[i],"*") != 0) a[i] = atof(cmdargs[i]);
  92.     }
  93. }
  94. else
  95.     printf("ip failed, there should be %d parameter values entered\n",ma);
  96. }
  97.  
  98. /* This function changes one of the parameters */
  99. void cp(command, inbuf, func, a, ma)
  100. char *command, *inbuf;
  101. int (*func)();
  102. double *a;
  103. int ma;
  104. {
  105. int i;
  106. char cmdargs[NUM_ARGS][30];
  107.  
  108. if(func == NULL)
  109.     printf("cp failed, you must select a function first\n");
  110. else if(parse(command,cmdargs) == 2 && atoi(cmdargs[0]) > -1
  111.         && atoi(cmdargs[0]) < ma)
  112.             a[atoi(cmdargs[0])] = atof(cmdargs[1]);
  113. else
  114.     help("cp",1);
  115. }
  116.  
  117. /* sp selects the parameters to vary */
  118. int sp(command, inbuf, func, lista, ma)
  119. char *command, *inbuf;
  120. int (*func)();
  121. int *lista, ma;
  122. {
  123. int i, mfit, max=0, min = 0;
  124. char cmdargs[NUM_ARGS][30];
  125.  
  126. if(func == NULL){
  127.     printf("sp failed, you must select a function first\n");
  128.     return 0;
  129. }
  130. else{
  131.     mfit = parse(command,cmdargs);
  132.     if(debug) printf("mfit: %d ma: %d\n", mfit, ma);
  133.     /* if mfit > ma, default to fitting all parameters */
  134.     if(ma < mfit){
  135.         mfit = ma;
  136.         for(i=0; i < ma; i++) lista[i] = i;
  137.         printf("sp failed, too many parameters selected\n");
  138.         if(debug)printf("mfit: %d ma: %d\n", mfit, ma);
  139.     }
  140.     else{
  141.         /* assign lista[i] to command arguments */
  142.         for( i = 0; i < mfit; i++) lista[i] = atoi(cmdargs[i]);
  143.         /* if any lista[i]'s are too big or too small, Default to all parameters */
  144.         for(i = 0; i < mfit; i++){
  145.             if( lista[i] > max) max = lista[i];
  146.             if( lista[i] < min) min = lista[i];
  147.         }
  148.         if(max > ma-1 || min < 0){
  149.             printf("sp failed, parameter out of range, all parameters fitted\n");
  150.             mfit = ma;
  151.             for(i=0; i < ma; i++) lista[i] = i;
  152.         }
  153.     }
  154.     }
  155. return mfit;
  156. }
  157.  
  158. int make_data(int (*func)(), int num_indep, double *a, int ma, char *command, char *inbuf){
  159.  
  160. char cmdargs[NUM_ARGS][30];
  161. char filename[80];
  162. FILE *stream;
  163. double *x, *xmin, *xmax, *xstep, y;
  164. int i,j, *fita;
  165. double  *dydx, *dyda;
  166. int *dydx_flag;
  167. int failed;
  168. double ydat = 0;
  169.  
  170. /* parse and check number of args entered on command line */
  171. if( (parse(command,cmdargs)) != 1 + 3*num_indep || func == NULL){
  172.     help("md",1);
  173.     return 1;
  174. }
  175.  
  176. if( (stream = fopen(cmdargs[0],"w")) == NULL){
  177.     printf("md error, cannot open file %s\n", cmdargs[0]);
  178.     return 1;
  179. }
  180.  
  181. if( num_indep > 2 ){
  182.     printf("md error, md not implemented for more then 2 independent parameters\n");
  183.     return 1;
  184. }
  185.  
  186. /* need to send all these to the function */
  187. dydx = dvector(num_indep);
  188. dydx_flag = ivector(num_indep);
  189. fita = ivector(ma);
  190. dyda = dvector(ma);
  191.  
  192. /* we only need function value, not derivatives */
  193. for(j = 0; j < num_indep; j++) dydx_flag[j] = 0;
  194. for(j = 0; j < ma; j++) fita[j] = 0;
  195.  
  196. x = dvector(num_indep);
  197. xmin = dvector(num_indep);
  198. xmax = dvector(num_indep);
  199. xstep = dvector(num_indep);
  200.  
  201. /* assign command line arguments */
  202. for( j = 0; j < num_indep; j++){
  203.     x[j] = xmin[j] = atof(cmdargs[3*j+1]);
  204.     xmax[j] = atof(cmdargs[3*j+2]);
  205.     xstep[j] = atof(cmdargs[3*j+3]);
  206.     if(debug)printf("xmin: %g xmax: %g xstep: %g\n", xmin[j], xmax[j], xstep[j]);
  207. }
  208.  
  209. if(num_indep == 1){
  210.     while( x[0] <= xmax[0] ){
  211.         failed = (*func)(x,a,&y,dyda,ma,0,fita,dydx_flag,dydx,ydat);
  212.         if(failed){
  213.             free(dydx);
  214.             free(dydx_flag);
  215.             free(fita);
  216.             free(dyda);
  217.             free(x);
  218.             free(xmin);
  219.             free(xmax);
  220.             free(xstep);
  221.         return 1;
  222.         }
  223.         fprintf(stream,"%g %g\n", x[0], y);
  224.         x[0] += xstep[0];
  225.     }
  226. }
  227.  
  228. if(num_indep == 2){
  229.     while( x[0] <= xmax[0] ){
  230.         x[1] = xmin[1];
  231.         while( x[1] <= xmax[1] ){
  232.             failed = (*func)(x,a,&y,dyda,ma,0,fita,dydx_flag,dydx,ydat);
  233.             if(failed){
  234.                 free(dydx);
  235.                 free(dydx_flag);
  236.                 free(fita);
  237.                 free(dyda);
  238.                 free(x);
  239.                 free(xmin);
  240.                 free(xmax);
  241.                 free(xstep);
  242.                 return 1;
  243.             }
  244.             fprintf(stream,"%g %g %g \n", x[0], x[1], y);
  245.             if(debug)printf("%g %g %g\n", x[0], x[1], y);
  246.             x[1] += xstep[1];
  247.         }
  248.         x[0] += xstep[0];
  249.     }
  250. }
  251.  
  252. fclose(stream);
  253. strcpy(inbuf,command);
  254.  
  255. free(dydx);
  256. free(dydx_flag);
  257. free(fita);
  258. free(dyda);
  259. free(x);
  260. free(xmin);
  261. free(xmax);
  262. free(xstep);
  263. return 0;
  264.  
  265. }
  266.  
  267. void est_errors(factor,func, data, order, num_indep, a, ndata, ma, chisqr)
  268. double factor;
  269. int (*func)();
  270. double **data;
  271. struct data_order order;
  272. int num_indep;
  273. double a[];
  274. int ndata, ma;
  275. double chisqr;
  276. {
  277. int i, j, failed,k;
  278. double tchisqr, tai, x1, x2, y1, y2, x;
  279. double perr, merr;
  280.  
  281. for(i = 0; i < ma; i++){  /* loop over parameters */
  282.     tai = a[i];              /* temporary a[i] */
  283.     tchisqr = chisqr;        /* temporary chisqr */
  284.     x = factor*chisqr;       /* chisqr we are looking for */
  285.  
  286.     /* change a[i] and use an iterative linear interpolation */
  287.     a[i] += 1e-7*tai;
  288.     failed = calc_yfit(func, data, order, num_indep, a, ndata, ma, &tchisqr);
  289.     x2 = chisqr;
  290.     x1 = tchisqr;
  291.     y1 = a[i];
  292.     y2 = tai;
  293.     a[i] = ( (x - x2)*y1 - (x - x1)*y2 )/(x1-x2);
  294.  
  295.     k = 0;
  296.     while( fabs(tchisqr - factor*chisqr) > 1e-5*fabs(chisqr) && k < 1000){
  297.     if(debug)printf("tai: %g a[%d]: %g chisqr: %g tchisqr: %g\n", tai, i, a[i], chisqr, tchisqr);
  298.     x2 = x1;
  299.     y2 = y1;
  300.     failed = calc_yfit(func, data, order, num_indep, a, ndata, ma, &tchisqr);
  301.     x1 = tchisqr;
  302.     y1 = a[i];
  303.     a[i] = ( (x - x2)*y1 - (x - x1)*y2 )/(x1-x2);
  304.     k++;
  305.     }
  306.     perr = a[i] - tai;
  307.     a[i] = tai;
  308.  
  309.     /* repeat the process for negative error, except begin with
  310.     assuming negative error is the same as the positive error */
  311.  
  312.     tchisqr = chisqr;
  313.     a[i] -= perr;
  314.     failed = calc_yfit(func, data, order, num_indep, a, ndata, ma, &tchisqr);
  315.     x2 = chisqr;
  316.     x1 = tchisqr;
  317.     y1 = a[i];
  318.     y2 = tai;
  319.     a[i] = ( (x - x2)*y1 - (x - x1)*y2 )/(x1-x2);
  320.  
  321.     k = 0;
  322.     while( fabs(tchisqr - factor*chisqr) > 1e-5*fabs(chisqr) && k < 1000){
  323.     if(debug)printf("tai: %g a[%d]: %g chisqr: %g tchisqr: %g\n", tai, i, a[i], chisqr, tchisqr);
  324.     x2 = x1;
  325.     y2 = y1;
  326.     failed = calc_yfit(func, data, order, num_indep, a, ndata, ma, &tchisqr);
  327.     x1 = tchisqr;
  328.     y1 = a[i];
  329.     a[i] = ( (x - x2)*y1 - (x - x1)*y2 )/(x1-x2);
  330.     k++;
  331.     }
  332.  
  333.     merr = -a[i] + tai;
  334.     a[i] = tai;
  335.  
  336.     printf("a[%d] = %g (+ %g, -%g)\n", i, a[i], fabs(perr), fabs(merr));
  337. }
  338.  
  339. }
  340.